library(Seurat)
library(Matrix)

readCDNARawMat <- function(prefix=NULL,data.dir) {
  if (!dir.exists(data.dir)) {
    stop(paste("Directory ",run," does not exist",sep=""))
  }
  if (!grepl("\\/$",data.dir))
  {
    data.dir <- paste(data.dir, "/", sep = "")
  }

  barcode.loc <- paste0(data.dir, "barcodes.tsv")
  gene.loc <- NULL
  gene.loc.1 <- paste0(data.dir,"genes.tsv")
  gene.loc.2 <- paste0(data.dir, "features.tsv")
  matrix.loc <- paste0(data.dir, "matrix.mtx")

  if (!file.exists(barcode.loc)) {
    stop("barcodes.tsv file is not existed!")
  }
  if (file.exists(gene.loc.1)) {
    print("find genes.tsv")
    gene.loc <- gene.loc.1
  }else if(file.exists(gene.loc.2)){
    print("find features.tsv")
    gene.loc <- gene.loc.2
  }else {
    stop("can't find genes.tsv or features.tsv")
  }

  if (!file.exists(matrix.loc)) {
    stop("matrix.mtx is not existed!")
  }

  data <- readMM(file = matrix.loc)
  cell.names <- readLines(barcode.loc)
  cell.names <- paste(prefix,cell.names,sep="")
  gene.names <- readLines(gene.loc)

  GeneName <- make.unique(names = as.character(unlist(lapply(gene.names,function(x)unlist(strsplit(x,"\t"))[2]))))
  rownames(data) <- GeneName
  colnames(data) <- cell.names
  return(data)
}

filterCellsNew <- function (x,cutoff=0.2) {
  #filter Dead Cells
  mito.genes <- grep(pattern = "^MT-", x = rownames(x = x), value = TRUE)
  if(length(mito.genes) == 0){
    mito.genes <- grep(pattern = "^mt-", x = rownames(x = x), value = TRUE)
  }
  nUMI <- Matrix::colSums(x)
  percent.mito <- Matrix::colSums(x[mito.genes, ])/nUMI
  deadcells <- which(percent.mito>cutoff)
  filter.cells <- x[,-deadcells]

  #filter Gene Exp is too big or too low
  x <- filter.cells
  num.gene <- colSums(x > 0)
  filterGene <- names(x=num.gene[which(x = (num.gene> 200 & num.gene < 2500))])
  filCellGeneMat <- x[,filterGene]

  x <- filCellGeneMat
  is.expr <- 0
  nUMI <- Matrix::colSums(x)
  nGene <- colSums(x > is.expr)
  lm.x <- lm(nGene~nUMI)
  lm.x.res <- scale(lm.x$residuals)
  lowComCells <- which(lm.x.res[,1]< -3)
  filCellGeneLow <- x[,-lowComCells]
  return(filCellGeneLow)
}
CDNAFMat <- readCDNARawMat(prefix="P_Tumor_", "cdna")
GCMat <- filterCellsNew(CDNAFMat, cutoff=0.2)
save(GCMat,file="GCMat.RData")



  
